# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python
source("preliminaries.R")
nbsimu <- 1000
alldf <- data.frame()
ps <- c(4, 8, 16, 32, 48, 64, 96, 128)

for (p in ps){
    load(paste("results-p", p, "T", p, "n", nbsimu, ".RData", sep = ""))
    df <- data.frame(cbind(basicestimvector, blockestimvectorNoPermut, blockestimvectorCircPermut, 
                           blockestimvectorRandPermut, blockestimvectorHalfRandPermut, blockestimvectorStratRandPermut))
    names(df) <- c("basic IMH", "same order", "circular", "random", "half-random", "stratified")
    df$p <- paste("p = ", p, sep = "")
    df$indexsimu <- 1:nbsimu
    alldf <- rbind(alldf, df)
}

alldf$p <- factor(alldf$p, levels = paste("p = ", ps, sep = ""))
test <- melt(alldf, id = c("indexsimu", "p"))
colourpalette <- hcl(h = seq(120, 0, length = 7) + 150, l = seq(60, 10, length = 7))

a <- ddply(test, .(p, variable), colwise(var, c("value")))

newa <- c()
for (index in 1:(dim(a)[1])){
    base <- subset(a, variable == "basic IMH" & p == a[index, 1])
    base <- as.numeric(base[3])
    newa <- c(newa, as.numeric((a[index, 3] - base) / base))
}
a$value <- newa
a <- subset(a, variable != "basic IMH")
a$variable <- factor(a$variable)
g <- ggplot(a, aes(x = variable, y = value))
g <- g + geom_bar(aes(size = 1, fill= variable))
g <- g + facet_wrap( ~ p, nrow = 2)
g <- g + xlab("") + ylab("") + opts(legend.position = "none")
g <- g + ylim(-0.5, 0)
g <- g + opts(axis.text.x = theme_text(angle = -45, size = 20, hjust = 0.2, vjust = 0.8))
g <- g + scale_fill_manual(values = colourpalette)
pdf(file = paste("comparisonPermutations-nbsim", nbsimu, ".pdf", sep = ""), width = 16, height = 5, useDingbats = FALSE)
print(g)
dev.off()

